{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.6 Error estimation & adaptive refinement\n",
    "\n",
    "\n",
    "In this tutorial, we apply a Zienkiewicz-Zhu type error estimator and run an adaptive loop with these steps:\n",
    "$$\n",
    "\\text{Solve}\\rightarrow\n",
    "\\text{Estimate}\\rightarrow\n",
    "\\text{Mark}\\rightarrow\n",
    "\\text{Refine}\\rightarrow\n",
    "\\text{Solve} \\rightarrow \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from ngsolve import *\n",
    "from netgen.geom2d import SplineGeometry\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometry\n",
    "\n",
    "The following geometry represents a heating chip embedded in another material that conducts away the heat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   point numbers 0, 1, ... 11\n",
    "#   sub-domain numbers (1), (2), (3)\n",
    "#  \n",
    "#\n",
    "#             7-------------6\n",
    "#             |             |\n",
    "#             |     (2)     |\n",
    "#             |             |\n",
    "#      3------4-------------5------2\n",
    "#      |                           |\n",
    "#      |             11            |\n",
    "#      |           /   \\           |\n",
    "#      |         10 (3) 9          |\n",
    "#      |           \\   /     (1)   |\n",
    "#      |             8             |\n",
    "#      |                           |\n",
    "#      0---------------------------1\n",
    "#\n",
    "\n",
    "def MakeGeometry():\n",
    "    geometry = SplineGeometry()\n",
    "    \n",
    "    # point coordinates ...\n",
    "    pnts = [ (0,0), (1,0), (1,0.6), (0,0.6), \\\n",
    "             (0.2,0.6), (0.8,0.6), (0.8,0.8), (0.2,0.8), \\\n",
    "             (0.5,0.15), (0.65,0.3), (0.5,0.45), (0.35,0.3) ]\n",
    "    pnums = [geometry.AppendPoint(*p) for p in pnts]\n",
    "    \n",
    "    # start-point, end-point, boundary-condition, domain on left side, domain on right side:\n",
    "    lines = [ (0,1,1,1,0), (1,2,2,1,0), (2,5,2,1,0), (5,4,2,1,2), (4,3,2,1,0), (3,0,2,1,0), \\\n",
    "              (5,6,2,2,0), (6,7,2,2,0), (7,4,2,2,0), \\\n",
    "              (8,9,2,3,1), (9,10,2,3,1), (10,11,2,3,1), (11,8,2,3,1) ]\n",
    "        \n",
    "    for p1,p2,bc,left,right in lines:\n",
    "        geometry.Append([\"line\", pnums[p1], pnums[p2]], bc=bc, leftdomain=left, rightdomain=right)\n",
    "\n",
    "    geometry.SetMaterial(1,\"base\")\n",
    "    geometry.SetMaterial(2,\"chip\")\n",
    "    geometry.SetMaterial(3,\"top\")    \n",
    "\n",
    "    return geometry\n",
    "\n",
    "mesh = Mesh(MakeGeometry().GenerateMesh(maxh=0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spaces & forms\n",
    "\n",
    "The problem is to find $u$ in $H_{0,D}^1$ satisfying \n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\lambda \\nabla u \\cdot \\nabla v = \\int_\\Omega f v \n",
    "$$\n",
    "\n",
    "for all $v$ in $H_{0,D}^1$. We expect the solution to have singularities due to the nonconvex re-enrant angles and discontinuities in $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=[1])\n",
    "u, v = fes.TnT()\n",
    "\n",
    "# one heat conductivity coefficient per sub-domain\n",
    "lam = CoefficientFunction([1, 1000, 10])\n",
    "a = BilinearForm(fes)\n",
    "a += SymbolicBFI(lam*grad(u)*grad(v))\n",
    "\n",
    "# heat-source in inner subdomain\n",
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(CoefficientFunction([0, 0, 1])*v)\n",
    "\n",
    "c = Preconditioner(a, type=\"multigrid\", inverse=\"sparsecholesky\")\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the linear system is not yet assembled above.\n",
    "\n",
    "### Solve \n",
    "\n",
    "Since we must solve multiple times, we define a function to solve the boundary value problem, where assembly, update, and solve occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP():\n",
    "    fes.Update()\n",
    "    gfu.Update()\n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    inv = CGSolver(a.mat, c.mat)\n",
    "    gfu.vec.data = inv * f.vec\n",
    "    Redraw (blocking=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "SolveBVP()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate\n",
    "\n",
    "We implement a gradient-recovery-type error estimator. For this, we need an H(div) space for flux recovery. We must compute the flux  of the computed solution and interpolate it into this H(div) space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "space_flux = HDiv(mesh, order=2)\n",
    "gf_flux = GridFunction(space_flux, \"flux\")\n",
    "\n",
    "flux = lam * grad(gfu)\n",
    "gf_flux.Set(flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = 1/lam*(flux-gf_flux)*(flux-gf_flux)\n",
    "Draw(err, mesh, 'error_representation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element-wise error estimator:** On each element $T$, set \n",
    "\n",
    "$$\n",
    "\\eta_T^2 = \\int_T \\frac{1}{\\lambda} \n",
    "|\\lambda \\nabla u_h - I_h(\\lambda \\nabla u_h) |^2\n",
    "$$\n",
    "\n",
    "where $u_h$ is the computed solution `gfu` and $I_h$ is the interpolation performed by `Set` in NGSolve.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 2.36461e-10\n",
      " 1.18698e-08\n",
      " 5.03613e-06\n",
      " 6.34682e-11\n",
      " 4.33101e-11\n",
      " 1.7484e-11\n",
      " 1.58609e-10\n",
      " 6.33751e-09\n",
      " 1.19293e-10\n",
      " 1.35243e-07\n",
      " 1.98437e-07\n",
      " 2.00923e-07\n",
      " 1.16473e-08\n",
      " 1.251e-07\n",
      " 2.2953e-07\n",
      " 5.73909e-09\n",
      " 3.385e-11\n",
      " 3.31297e-10\n",
      " 3.32684e-10\n",
      " 3.89806e-07\n",
      " 2.26934e-07\n",
      " 2.60348e-07\n",
      " 3.3506e-07\n",
      " 8.52564e-10\n",
      " 2.40196e-06\n",
      " 3.57631e-06\n",
      " 2.45321e-08\n",
      " 1.29118e-10\n",
      " 8.8645e-08\n",
      " 3.87456e-08\n",
      " 1.6248e-07\n",
      " 7.55763e-08\n",
      " 3.12331e-06\n",
      " 3.79783e-10\n",
      " 1.90566e-07\n",
      " 1.04137e-09\n",
      " 6.8768e-07\n",
      " 7.2601e-07\n",
      " 3.8487e-07\n",
      " 1.43514e-10\n",
      " 1.35152e-10\n",
      " 6.55815e-07\n",
      " 1.82416e-07\n",
      " 7.13649e-08\n",
      " 9.56596e-09\n",
      " 6.74125e-09\n",
      " 3.28339e-10\n",
      " 2.90424e-09\n",
      " 5.48405e-08\n",
      " 1.15509e-08\n",
      " 1.45004e-08\n",
      " 2.01172e-08\n",
      " 1.37504e-08\n",
      " 1.13422e-07\n",
      " 3.51327e-08\n",
      " 1.27836e-10\n",
      " 5.86447e-12\n",
      " 4.22009e-12\n",
      " 3.1807e-10\n",
      " 5.72781e-13\n",
      " 1.2421e-11\n",
      " 5.4183e-12\n",
      " 2.41244e-10\n",
      " 6.65466e-07\n",
      " 8.8388e-08\n",
      "\n"
     ]
    }
   ],
   "source": [
    "eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "print(eta2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above values, one per element, lead us to identify elements which might have large error.\n",
    "\n",
    "\n",
    "### Mark \n",
    "\n",
    "We mark elements with large error estimator for refinement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maxerr =  5.036125830814317e-06\n"
     ]
    }
   ],
   "source": [
    "maxerr = max(eta2)\n",
    "print (\"maxerr = \", maxerr)\n",
    "\n",
    "for el in mesh.Elements():\n",
    "    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)\n",
    "    \n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refine & solve again \n",
    "\n",
    "Refine marked elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.Refine()\n",
    "SolveBVP()\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automate the above steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []    # l = list of estimated total error\n",
    "\n",
    "def CalcError():\n",
    "\n",
    "    # compute the flux:\n",
    "    space_flux.Update()      \n",
    "    gf_flux.Update()\n",
    "    flux = lam * grad(gfu)        \n",
    "    gf_flux.Set(flux) \n",
    "    \n",
    "    # compute estimator:\n",
    "    err = 1/lam*(flux-gf_flux)*(flux-gf_flux)\n",
    "    eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "    maxerr = max(eta2)\n",
    "    l.append ((fes.ndof, sqrt(sum(eta2))))\n",
    "    print(\"ndof =\", fes.ndof, \" maxerr =\", maxerr)\n",
    "    \n",
    "    # mark for refinement:\n",
    "    for el in mesh.Elements():\n",
    "        mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 445  maxerr = 2.511445609726963e-06\n"
     ]
    }
   ],
   "source": [
    "CalcError()\n",
    "mesh.Refine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the adaptive loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndof = 664  maxerr = 9.300425565031251e-07\n",
      "ndof = 1192  maxerr = 3.414646836906428e-07\n",
      "ndof = 1834  maxerr = 2.1472585059819264e-07\n",
      "ndof = 2167  maxerr = 9.835056165137945e-08\n",
      "ndof = 2839  maxerr = 4.917040165457557e-08\n",
      "ndof = 3466  maxerr = 2.458740959858597e-08\n",
      "ndof = 3805  maxerr = 1.2294699792242802e-08\n",
      "ndof = 4450  maxerr = 6.145986051434467e-09\n",
      "ndof = 5257  maxerr = 3.071914052546866e-09\n",
      "ndof = 5713  maxerr = 1.5353408223467698e-09\n",
      "ndof = 6475  maxerr = 7.673155137357938e-10\n",
      "ndof = 7165  maxerr = 3.8347357373724957e-10\n",
      "ndof = 7900  maxerr = 1.9164536414781034e-10\n",
      "ndof = 8719  maxerr = 9.577670728615699e-11\n",
      "ndof = 9853  maxerr = 4.7867485757047704e-11\n",
      "ndof = 10810  maxerr = 2.392324427586057e-11\n",
      "ndof = 11884  maxerr = 1.1956352408132547e-11\n",
      "ndof = 13513  maxerr = 5.975487707207399e-12\n",
      "ndof = 15553  maxerr = 2.986452644572587e-12\n",
      "ndof = 17776  maxerr = 1.4925627796160417e-12\n",
      "ndof = 20404  maxerr = 7.459783254409289e-13\n",
      "ndof = 23653  maxerr = 3.7281573519588415e-13\n",
      "ndof = 27394  maxerr = 1.863284713063423e-13\n",
      "ndof = 31690  maxerr = 9.312791086715527e-14\n",
      "ndof = 36811  maxerr = 4.6543377383552036e-14\n",
      "ndof = 42022  maxerr = 2.3263412632625325e-14\n",
      "ndof = 49306  maxerr = 1.1626515614277681e-14\n",
      "ndof = 57589  maxerr = 5.81163901405017e-15\n"
     ]
    }
   ],
   "source": [
    "while fes.ndof < 50000:  \n",
    "    SolveBVP()\n",
    "    CalcError()\n",
    "    mesh.Refine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot history of adaptive convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl81PW1//HXyc6SBIGETQIoCKIsIosotoK1V6vg0l6toFZbd7G/ri6trW2vt/Z6W229aJW2ahX3pQoqaluhgFAMi4jsIBBAMIQl7CHJnN8fk2DEQCbJ7PN+Ph7zIPPNdyZHHXP4fs/5nI+5OyIiIqFKi3UAIiKSWJQ4RESkUZQ4RESkUZQ4RESkUZQ4RESkUZQ4RESkUZQ4RESkUZQ4RESkUZQ4RESkUZQ4RESkUTJiHUA4mdloYHRubu51J5xwQqzDERFJKPPnzy9z94KGzrNknFU1ePBgnzdvXqzDEBFJKGY2390HN3SeblWJiEijKHGIiEijKHGIiEijKHGIiEijJFXiMLPRZjaxvLy8Sa8v3XWASx+dQ+nuA2GOTEQkeSRV4nD3Ke5+fX5+fpNe/+A/V1G8bjsP/mNVmCMTEUkeSbWOo6l63zWViqrAoeeT5pYwaW4J2RlprLjnvBhGJiISf5LqiqOpZt42kjEDO5Nmnx3r0qYFP/laH7btqYhdYCIicUhXHEBhXg652Rk4kJluVFY7ZXsquHvyUn4xZSkDu7ZhZO9CRvYu5KTOeaTVzTAiIilGiaNG2Z4Kxg3rxtihRTzzfglbdx1g/KhevLu8lGkrSnngHyu5/+8rKcjN5qwTChjZp5ARvdqTl5MZ69BFRKJKI0dCtG1PBf9auZVpK7YyY+VWyvdXkpFmDO5+TPBqpE8hvQpbY6arERFJTKGOHEmqxFE75LBnz57XrVoVuc6oquoAH2zYWXM1spVlm3cBwbrIyD4FjOxdyPDj29EySxd0IpI4UjJx1Ir2kMPN5fuZvmIr05aXMmt1GfsOVpOVkcZpx7VjVO/gba1u7VoBwbUi459dyISxp1CYmxO1GEVEGqLEEaPpuBVV1RSv3cG0FcHayMdb9wJwXPtWjOxTyLqyvby7opRxQ4u45+J+MYlRRKQ+ShxxMlZ9/ba9TFteyq9eX0qgnn/VmenGgp+dQ66K7CISY0occZI4apXuOsAvpyzl70u3cLDaMaD233yawYmd8hjSvS2Dux/DkO5t6ZCn21giEl2hJg5Vb6OkMC+HNi0zqQw42RlpHKwOcNmpXRk9oDPF67Yzb/12ni/ewBOz1wFQ1LbloSQypPsxHF+gji0RiQ9KHFH0hbUiuw8wold7RvRqD0BldYCln+wKJpJ1O5ixciuvLNgEwDEtMxlck0QGd2/LyZ3zycrQwn8RiT7dqopj7s7asr3MW7ej5qpkB2vLgsX2nMw0BnZtU3N7qy2DitqoTiIizaIaRxIkjvqU7j7A/HU7KF63g3nrt7Pkk11UB1x1EhFpNiWOJE0ch9tTUcUHJTsP1UkWrN/J/spq4PA6SVuOL2hVb51Ea0tEBFQcTxmtszOOWif514r66yRDurflpJo6Sd19SLS2REQaklRXHNEaOZJIDq+TFK/bzrpt+476Gu1DIpKadKsqRW5VNUVtneRfK7cy9aMtlO+vBMAMBhW14cFvDqLLMS1iHKWIRFuoiUP9nCmoMDeH8/p14jdf788F/TthBhlphjvMX7+TCx96j9+9s4LN5ftjHaqIxCHVOFLc59aWzF3Psi27adMikwnTVvPw9DV8tW8HrhzejeHHtdMCRBEBdKtKjmDD9n1Mmrue54s3sHNfJb0KW3PV8G5cPOhYWmfr7xsiyUg1DiWOsDhQWc2URZ/w5Jz1LN5UTuvsDC4Z1IWrhnejZ2FurMMTkTBS4lDiCCt354MNO3lqznpe/3AzB6sDnH58O64a3o2vnNiBjHSVy0QSnRKHEkfEbNtTwXPFG3hmbgmbdu6nU34O44YVcdmQIgpys2Mdnog0kRKHEkfEVQecfy77lCfnrGfW6jIy042v9evEVcO7M6iojYrpIglGK8cl4tLTjK+e1JGvntSRNVv38NSc9bw8fyOvffAJJ3XO46rh3RgzoAststJjHaqIhJGuOCSs9lZU8beFm3hqznpWfLqb/BaZXDr4WK44rRvd2rXSXCyROKZbVUocMeXuvL92O0/OWc9bS7YQcOfLJxRgwPSVW7XnukgcUuJQ4ogbn+46wOm/eZfqejZd11wskfihkSMSNzrk5TDnjlGM7t+JzPTPCub9j81nxo9HxjAyEWmKuE8cZnaimT1iZi+Z2U2xjkeapjAvh7wWmVQF/NCWtx9uLOeHLy5i007NxBJJJBFNHGb2mJmVmtlHhx0/18xWmNlqM7vjaO/h7svc/UbgUuCMSMYrkVU7F+vVm8/gimFF9O2Ux4KSHZz7wAxemLeBZLxtKpKMIlrjMLMvAXuAJ9395Jpj6cBK4BxgI1AMXA6kA/ce9hbfdvdSMxsD3AQ85e7PNPRzVeNIHCXb9vHjlxYxd+12zu5TyL2X9KNQW96KxERc1DjcfQaw/bDDQ4HV7v6xux8EngMudPfF7n7BYY/SmveZ7O7nAeMiGa9EX1G7ljx73Wn8/IK+zFpdxjkPzGDyok/4tHw/lz46h9LdB2IdoogcJhYLALsAG+o83wgMO9LJZnYWcAmQDbx5lPOuB64HKCoqCkecEiVpaca3R/Tgy70L+NGLi/juswspatuCDTv2aztbkTgU9yvH3X06MD2E8yYCEyF4qyqyUUkkHF/QmqWf7AKgZHuwYD5pbgmT5paobVckjsSiq2oT0LXO82NrjjWbmY02s4nl5eXheDuJgZm3jWTMwM5kZ3z20WzfKotJ1x7xolREoiwWiaMY6GVmPcwsC/gmMDkcb+zuU9z9+vz8/HC8ncRAYV4OudkZHKwOkJ2RhgHlByoZ96e53P/3lRyorI51iCIpL6K3qszsWeAsoL2ZbQTudve/mNl44G2CnVSPufuSSMYhieVz29m+X8LGHfvIb5HJg/9cxeQPNnHPRf0Y0at9rMMUSVlJNXLEzEYDo3v27HndqlWrYh2OhNmsVWXc9epi1m3bx4UDO3PX+X21/4dIGGlWldZxJKUDldU8PH0Nj0xfQ3ZmGref24exQ4tIS9PeHyLNFRfrOETCLScznR+ccwJTv3cmJ3fO565XP+Lrj8w+1I0lIpGnxCEJ6fiC1jxz3TDuv3QAJdv2MXrCLP77jaXsraiKdWgiSS+pEofacVOLmXHJoGP55w+/zKWDj+VPM9dyzv3/4p0lWwAo3XVAq89FIkA1Dkka89Zt56d/+4gVn+7mnL4daJ2dzqsffKJNo0RCpOK4EkdKqqwO0Odnb2nTKJEmUHFcUlJmehpz7hjFV04spLbRKs3ggv6dmHm7No0SCYekShyqcQgEV593yMvBgfQ0I+Awc1UZ+yq06lwkHJIqcWjkiNSqXX0+ZfwIzu5TyN6KKsZMmMW05aWxDk0k4anGISlhw/Z93PDUfJZt2cX3v3IC40f21KJBkcOoxiFSR9e2LXn5ptO5aGAX7v/7Sq5/aj67DlTGOiyRhKTEISmjRVY69186gF+M7sv0FaVcNOE9Vn26O9ZhiSScpEocKo5LQ8yMq8/owdPXDmPXgUoueug9pi7eHOuwRBJKUiUOFcclVMOOa8frt55Jrw653PT0Au57aznVAddqc5EQxP3WsSKR0jE/h+dvOI1fTF7Kw9PX8NEnu+iQm03xuu3a61zkKNRVJQL0/MmbVGm1uaQ4dVWJNMLsO0ZxZp1dBTPSjDEDtNpcpD4hJQ4zG2Fm19R8XWBmPSIbVtOoOC5NVZiXQ1HblpgFR5RUBZz3Vm9jz4Eq1T1EDtNg4jCzu4HbgTtrDmUCkyIZVFOpOC7NUbva/PVbz+TMXu0o31/J1x6cyS3PLDxU9xCREGocZvYBcAqwwN1PqTn2obv3j0J8TaIah4TDCT+dysHqwBeOq+4hySqcNY6DHswuXvPGrZobnEgimHX7SMYM6ExmenA0iQH/cVIH1T0k5YWSOF4ws0eBNmZ2HfAP4M+RDUsk9grzcsjNyaAq4GSmGw7MWlXGXk3ZlRTXYOJw998CLwEvA72Bn7v7g5EOTCQe1NY9XrtlBOed3JGD1QG+8cfZfLRJDRiSukKpcfyPu9/e0LF4ohqHRMrHW/dw5V/ep3x/JROvOpXTj2/f8ItEEkQ4axzn1HNMlUFJSccVtOblm06nU34OVz9WzFsfac6VpJ4jJg4zu8nMFgO9zezDOo+1wIfRCzF0Wsch0dAxP4cXbxzOSV3yuPnpBTz7fkmsQxKJqqNdcTwDjAYm1/xZ+zjV3a+IQmyNpnUcEi1tWmbx9LXDOLNXAXe+spiHpq3m0/L9WigoKeGIicPdy919nbtf7u7rgf0EW3Jbm1lR1CIUiVMtszL401WDuXBgZ/737RV86/FiLRSUlNDgdFwzGw3cD3QGSoFuwDLgpMiGJhL/sjLSeOujLQAs3xLcFGrS3BImzS3RQkFJWqEUx+8BTgNWunsP4Gzg3xGNSiSBzLwtuFAwo2YP83QzLhzYWQsFJWmFkjgq3X0bkGZmae4+DWiwXUskVdQuFKx2Jz3NqHZnxZbdFLTOjnVoIhERSuLYaWatgRnA02b2B2BvZMMSSSy1CwUnjz+D3h1as3zLbh5QrUOSVCgLAFsBBwiO6hkH5ANP11yFxCUtAJRYCgScO175kBfmbeQnX+vD9V86PtYhiYQk1AWADRbH3X1vzRvmAVPCEJtIUktLM+69pD97K6r59ZvLaZ2dydhhakSU5BFKV9UNwC8JXnUECF55OHBcZEMTSVzpacYDlw1kf2U1P311MS2z0rnolC6xDkskLEKpcfwIONndu7v7ce7ew93jMmlo5bjEk6yMNB4eN4hhPdrywxcX8c6SLbEOSSQsQkkca4B9kQ4kHLRyXOJNTmY6f/7WEPp1yWf8MwuZtaos1iGJNFsoieNOYLaZPWpmD9Y+Ih2YSLJonZ3BE9cM4biCVlz35Dzmr9+ufcwloYWSOB4F3iW46G9+nYeIhKhNyyye+s4wOubncPXjxdw9eYnGk0jCarA4DmS6+w8iHolIkivIzWbTjv0crA4wtWZMicaTSCIK5Ypjqpldb2adzKxt7SPikYkkoVm3j+TsPoWHnudkpGk8iSScUK44Lq/58846x9SOK9IEhXk5dMzPOdTTfqAqQOvsDApzc2IdmkjIQlkA2CMagYikirI9FYw7rRv5LTJ4aNoaZq9Rp5UkliMmDjMb5e7vmtkl9X3f3V+JXFgiyevRKz+b6LB9byXPvl/ClEWfMHpA5xhGJRK6o11xfJlgN9Xoer7ngBKHSDP9csxJrPp0Nz9+aRFtWmbyf++uZsLYU3TrSuJaKEMOe7j72oaOxRMNOZREsnV3BWMmzGLX/kr2VVYzbmgR91zcL9ZhSQoKdchhKF1VL9dz7KXGhyQi9RnxP++yufwAew9W4x5s0e1+xxv0vmtqrEMTqdfRahx9CG4Pm39YnSMP0HW0SJjMvG0k97y5jLcWb+FgdQCA/sfm8+ertF+axKejXXH0Bi4A2hCsc9Q+BgHXRT40kdRQmJdDbnYGlYEAWRnB/yU/3FjO3ZOXUL6/MsbRiXzREROHu7/m7tcAF7j7NXUe33X32VGMETNrZWbzzOyCaP5ckWip3UHw1ZvP4IphRZzQoTV/X/opX/vDTBaU7ADQfCuJG6EUx+8D7gH2A28B/YHvu/ukBt/c7DGCVy2l7n5ynePnAn8A0oE/u/tvGnifXwF7gKXu/npDP1fFcUkGC0p28N1nF7K5/AA/+mpvNu7YxzPvl6h4LhETanE8lMTxgbsPNLOLCSaBHwAz3H1ACEF8ieAv/CdrE4eZpQMrgXOAjUAxwdXp6cC9h73Ft4EBQDuCdZUyJQ5JJeX7KznlV+8QqOd/U823knAL29axQGbNn+cDL7p7uZmFFIS7zzCz7ocdHgqsdvePawJ9DrjQ3e8lmJg+x8zOAloBfYH9ZvamuwdCCkAkweW3yGTOHaO4/qn5LNoY3KAsOyONc0/uyE/PPzHG0UmqCiVxTDGz5QRvVd1kZgUEt5Ftqi7AhjrPNwLDjnSyu/8UwMyuJnjFUW/SMLPrgesBioq0v7Mkjw75LTi5Sz4fbizHgYqqAHsOVGmRoMRMg+s43P0O4HRgsLtXEtwN8MJIB1ZPHE8c7TaVu09098HuPrigoCCaoYlEXO18q8evGUJ+iwzeXV7KC8UbGn6hSAQ0mDjMrCVwM/DHmkOdgeY0mG8CutZ5fmzNMRE5gkevHMw9F53MyN6FzLx9FCN6tee2lz/kvreWs2XnfnVbSVSFsnL8ceAgwasOCP6Sv6cZP7MY6GVmPcwsC/gmMLkZ73eImY02s4nl5eXheDuRuJSXk8ljVw/h8qFFPDx9DZdO/Ld2E5SoCiVxHO/u9wGVAO6+DwipOm5mzwJzgN5mttHMvuPuVcB44G1gGfCCuy9pUvSHcfcp7n59fn5+ON5OJG5lpqfxyoKNAJRs36dRJRJVoRTHD5pZC4ITcTGz44GKUN7c3S8/wvE3gTdDDVJEvqh2VMnUxZuprA72647qU8hvvq41HhJZoVxx3E1w4V9XM3sa+CdwW0SjaiLdqpJUUjuqpCrgZKYHbwLMXLmV5Zt3xzgySXahdFX9HbgEuBp4lmB31fTIhtU0ulUlqaZ2VMlrt4zgklO6kJ2ZzjVPFDPp3+s1okQipsGV45872ewX7v6LyIUTHlo5LqlqT0UVtz6zgGkrttKnYy4rPt2tESUSsnDux1HXmCbGIyJR0Do7g9lrtgGwfMtuFc0lIhqbOEKbNRIjqnGIBIvmYwZ2JqOm7pFmMLp/J2bePjLGkUmyOGriMLN0M/t+nUOnRjieZlGNQ+Szonl1wMlIMwIenLSb3yKz4ReLhOCoicPdqwlOrq19ruGCIgmgtmg+efwIhnZvy6adB7jl6YUcrNL/wtJ8oYxVf4DghNzngb21x919QWRDazoVx0U+76k56/jZa0s4p28H7h7dlx+8sIgJY0/RoET5nHCOVR9Y8+ev6hxzYFRTAoskMxsNjO7Zs2esQxGJK1cO744DP39tCcs272LTzv08+I9V6raSJmlUO26i0BWHyBf1vmsqFfXcqtKGUFIrbO24ZpZvZvfX7Pk9z8x+Z2aqPoskmNpuq8w63Vbn91O3lTReKO24jwG7gUtrHrsITswVkQRSd0RJbbfVnI/LSAtxR0+RWqFOx73b3T+uefwSOC7SgTWF1nGIHF3dbqtRvQvYua+Sb/xxNhu274t1aJJAQumqmgP82N1n1Tw/A/ituw+PQnxNohqHSGjmr9/Bt58oJisjjSe/PZQTO+XFOiSJoXCOHLkReMjM1pnZOmACcEMz4xOROHBqt2N46cbhZKQZlz46h7kfb9NwRGlQQyvH04De7j4A6A/0d/dT3P3DqEQnIhHXq0MuL990OoW52Vz52Pv86KVF2lFQjiqUW1XzQrl0iSe6VSXSeCfcNbXeleVq100d4bxV9Q8z+5GZdTWztrWPMMQoInFk1m0jOb9fR9Jqmqyy0o0LB3ZWu658QSgrxy+r+fOWOsecOOys0spxkaYrzMuhTcssnOAY7IPVzv6D1RpLIl8QSo3jCnfvcdgj7pIGaDquSHPVtus+c90wcrMzeHd5KatLtRWtfF4oNY6F7n5KlOIJC9U4RJpvw/Z9XPzwbLIz0njl5tPpkKcrj2QXzhrHP83s62ZaXiqSSrq2bckT1wxh576DfOux91mzdY/adAUILXHcALwIHDSzXWa228x2RTguEYkDJ3fJ55ErT2V16R7GTvy32nQFCKE47u650QhEROLTtX+dR1XA+XR3BRDcw3zS3BK16aawUKbjmpldYWY/q3ne1cyGRj40EYkHh0/VBTizV3u16aawUG5VPQwMB8bWPN8DPBSxiEQkrtSdqlubPGatKuONDzfj7hpRkoJCWccxzN0HmdlCAHffYWZZEY6rSbSOQyQyatt0xw4t4onZa5m2vJRfTlnK/PU7aJWVcaj2oR0FU0Mo7bhzgdOB4poEUgC8E88tumrHFYmsQMDpdddUqgNf/P2h2kfiCmc77oPA34BCM/tvYBbw62bGJyIJLC3NmHPHKE4/vt2hYxpRkjpC6ap62szmA2cTnERwkbsvi3hkIhLXCvNy6NG+FXPWbAOCI0rWlu2lfavsGEcmkRbKFQfuvtzdH3L3CUoaIlKrbE8F407rxt9uPp0e7Vrx4cZybn56AXsqqmIdmkRQgzWORKQah0j0uTt/mbWWX7+5jJ6FrZl45WC6t28V67CkEcJZ4xARaZCZce2Zx/Hkt4dRuruCMRNmMWPlVrXrJiElDhEJqxG92jP5lhF0btOCqx9/nxsnzdeokiTTpMRhZovDHYiIJI+idi1ZW7aXgMOCkp24B0eVdL/jDXrfNTXW4UkzHbGryswuOdK3gI6RCUdEksXM20ZyzxvLeHPxZqpq1nuc2as9v7t0QIwjk+Y6Wjvu88DTBHf7O1xcDubXynGR+FGYl0NuTgbVHhxVUlntzFxVxvPvb+DmkT1JT9NODYnqiF1VNWs3vuXuH9XzvQ3u3jXSwTWVuqpE4sMNT82jIDeHsUOL+OvstfxrVRlbyg8w/Lh2PHDZQDrmx+XfQVNWqF1VR0scZwLr3b2knu8Ndve4/c2sxCESn9ydl+Zv5OevLSEnM43f/ucAzj6xA6W7DjD+2YVMGHuK9jiPoWYnjkSmxCES39Zs3cOtzyxk6eZdXH16dyqqqnmueAPjhhZpUGIMheOK4/+ov74BgLt/t+nhRZYSh0j8q6iqpu/P3qa6nt9BGpQYG+FYADgPmF/zGFPn69qHiEiTZWekM+fOUQzr0fbQMQ1KTAxH7Kpy97/Wfm1m36v7XEQkHArzcuhZ2Jr3124HgoMSP966R4MS41yoCwCTrxAiInGh7qDE49q3YvGmXXznr8WU76uMdWhyBCEVx81sgbsPikI8YaEah0hicneenlvCL6csoVN+Cx654lT6ds6LdVgpo9k1DjPbbWa7zGwX0L/269rjYY1WRITgoMQrTuvG8zcMp6Kqmkv++B6vLtykQYlx5oiJw91z3T2v5pFR5+tcd9dfAUQkYgYVHcPrt55J/2Pb8L3nP+DKv8zVoMQ40uAOgCIisVCQm82iDTsBWPHpHiA4KHHS3BK168aYxqqLSNyaedtIxgzsTGb6Z3Ot+nTMZdqPzopdUKLEISLxqzAvh9zsDKoCTlZG8NfV8i27ufav8/hoU3mMo0tdcZ84zOwsM5tpZo+Y2VmxjkdEoqtsTwXjhnXj1ZvP4IrTujGwaxu27qngwofe47dvr6CiqhpABfQoimiNw8weAy4ASt395DrHzwX+AKQDf3b33xzlbRzYQ3CU+8YIhisicejRKz/rDr3nouCvkZ37DvJfry9jwrTVvL1kC/d9oz8vz994qICueVeRFdEhh2b2JYK/9J+sTRxmlg6sBM4hmAiKgcsJJpF7D3uLbwNl7h4wsw7A/e4+rqGfq3UcIqlh2opSrnm8uN7vqYDeeOGYVdVs7j4D2H7Y4aHAanf/2N0PAs8BF7r7Yne/4LBHqbsHal63A9AcAhE5ZGTvQv75wy/TrW3LQ8eyM9I07yrCYlHj6AJsqPN8Y82xepnZJWb2KPAUMOEo511vZvPMbN7WrVvDFqyIxLfjC1ozold7avuuKqoC7Nh7UPt6RFDcF8fd/RV3v8HdL3P36Uc5b6K7D3b3wQUFBVGMUERirXbe1eNXD6FdqyxmrCrjV1OWcrAq0PCLpdFisQBwE1B329lja46JiDRJ3QL67DtHce+by3nsvbUs3LCDCWMH0aVNixhGl3xiccVRDPQysx5mlgV8E5gcjjc2s9FmNrG8XP3dIqkqOyOdX4w5iYfGDmLVp3s4/8GZTFteCqhlN1wimjjM7FlgDtDbzDaa2XfcvQoYD7wNLANecPcl4fh57j7F3a/Pz88Px9uJSAI7v38nJo8/g455OVzzRDH3vbWc3/9jpWZehYH2HBeRpHagspq+P3+LQD2/6tSy+3lx0Y4bbbpVJSKHy8lM5993ns0pXdscOpaVrpbd5kiqxKFbVSJSn8K8HPp2zsMAAw5WB1hXtldb1DZRUiUOEZEjqW3Zfemm4XRv15JFG8u55olitu89GOvQEk5S1TjMbDQwumfPntetWqXil4jUr3aL2l9NWUr71llMGDeIQUXHxDqsmEvJGoduVYlIKGq3qH35ptNJTzcufWQOf5m1FndXy24IkipxiIg0Rr9j83l9/JmM7FPIf72+lJufXsBv31mhlt0GJNWtqlpqxxWRxnB3ev50KtX19OymUstuSt6qUjuuiDSFmTHnjlGM6Nn+0LHMdFPL7hEkVeJQjUNEmqowL4du7VpiBmZQWe0s3lhOXk5mrEOLO0mVOEREmqN2m9op40fQr0seH5ft5eKHZ7OubG+sQ4srqnGIiBzBtBWlfP/5D6iudv73Pwdw7skdYx1SRKVkjUNEJJxG9i7k9VtHcFxBK26cNJ9fv7mMyupAyrfsJlXiUHFcRMLt2GNa8sKNw7nytG5MnPExY//0b37z1vKUbtnVrSoRkRD1/MmbVCVxy65uVYmIhNnsO0YxqnfBof3NM9KM0f07pVzLrhKHiEiICvNy6NSmBRikGVQFnPfWbGNvRXWsQ4sqJQ4RkUaobdl9/dYz+VKv9uzaX8n5D87kheINJOOt//qoxiEi0gyby/fzg+cXMefjbZx3ckfuvaQfbVpmxTqsJknJGoe6qkQk2jrlt2DStcO447w+/H3pp5z7+5nMXlOW1C27uuIQEQmTxRvL+X/PLWTttr2c2DGPZVt2MW5oEfdc3C/WoYUk1CuOjGgEIyKSCvodm8+mnftxh6WbdwEwaW4Jk+aWJE3LLiTZrSoRkVibedtIxgzsTGa6HTrWr0seM36cPC27ShwiImFUmJdDbnYGVQEnKyP4K3bxpl187/kP2LhjX4yjCw8lDhGRMKtt2X315jO4YlgRJ3XO48ONOzn39zN5cV7it+2qOC4iEgVlXpbkAAAFxklEQVQbtu/jhy8u4v212zmnbwfuvaQfgYAz/tmFTBh7CoW5ObEOMeTieFIlDjMbDYzu2bPndatWpebwMRGJX4GA89h7a7nv7RW0zs7g5M55zFxdFjedVymZOGrpikNE4tkJP53KwerAF47HuvMqJRcAiogkglm3j2R0/05kpH3WeXX68e0SZliiEoeISJQV5uWQ1yKTavdDbbuz12zjj9PXcKAy/gcmKnGIiMRAbefVa7eM4PIhXSlq24LH31vH+Q/OZNGGnbEO76hU4xARiROzVpXx45cWUbq7gltH9eSWkT3JTI/e3+9V4xARSTAjerXnre99iTEDOvP7f6zi63+czerSPXE3MFGJQ0QkjuS3yOSBywby8LhBlGzfx/kPzuTmpxfE1R7nulUlIhKnot22q1tVIiIJbtbtIxkz4PCBifkxH5ioxCEiEqcK83LIzakZmJheOzCxnFufW8j6bXtjFldSJQ7tACgiyebQwMRbPhuYuGzzLv7j9zN4bNZaAoHolxtU4xARSTCby/fzk1cWM23FVoZ0P4b7vjGAVlnpzR6YqBqHiEiS6pTfgseuHsLv/nMAK7bs5tzfz+CmSfOj1nmlKw4RkQQWzs4rXXGIiKSAwzuvsjPSuHBg54gOTFTiEBFJYHU7r7Iz0jhYHSA3OyOiG0NlROydRUQkKmo7r8YOLeKZ90vYGuHRJKpxiIgIoBqHiIhEiBKHiIg0ihKHiIg0ihKHiIg0ihKHiIg0ihKHiIg0SlK245rZVmD9UU7JB5oyQrexrwvl/Oaec7TvtQfKGnjveNHU/yax+jnR+AyFem5D5+nzE38/J14/P93cvaDBV7h7yj2AidF4XSjnN/ecBr43L9b/riP93yRWPycan6FQz23oPH1+4u/nJNLnp75Hqt6qmhKl14VyfnPPaeo/S7yJ1j9HuH5OND5DoZ7b0Hn6/MTfz0mkz88XJOWtKgkys3kewipQkfro8yNHkqpXHKliYqwDkISmz4/US1ccIiLSKLriEBGRRlHiEBGRRlHiEBGRRlHiSBFmdqKZPWJmL5nZTbGORxKPmbUys3lmdkGsY5HYUuJIYGb2mJmVmtlHhx0/18xWmNlqM7sDwN2XufuNwKXAGbGIV+JLYz4/NW4HXohulBKPlDgS2xPAuXUPmFk68BBwHtAXuNzM+tZ8bwzwBvBmdMOUOPUEIX5+zOwcYClQGu0gJf5oz/EE5u4zzKz7YYeHAqvd/WMAM3sOuBBY6u6Tgclm9gbwTDRjlfjTyM9Pa6AVwWSy38zedPdAFMOVOKLEkXy6ABvqPN8IDDOzs4BLgGx0xSFHVu/nx93HA5jZ1UCZkkZqU+JIEe4+HZge4zAkwbn7E7GOQWJPNY7kswnoWuf5sTXHREKhz480SIkj+RQDvcysh5llAd8EJsc4Jkkc+vxIg5Q4EpiZPQvMAXqb2UYz+467VwHjgbeBZcAL7r4klnFKfNLnR5pKQw5FRKRRdMUhIiKNosQhIiKNosQhIiKNosQhIiKNosQhIiKNosQhIiKNosQhEiVm1v3wEeZHOO9ZM/vQzL4fjbhEGkuzqkTiiJl1BIa4e89YxyJyJLriEAmDmquJZWb2JzNbYmbvmFkLMzvVzBaZ2SLgljrn55jZ42a22MwWmtnImm+9A3Qxsw/M7MyY/MOINECJQyR8egEPuftJwE7g68DjwK3uPuCwc28B3N37AZcDfzWzHGAMsMbdB7r7zCjGLhIyJQ6R8Fnr7h/UfD0f6A60cfcZNceeqnPuCGASgLsvB9YDJ0QpTpFmUeIQCZ+KOl9XA+1jFYhIJClxiETOTmCnmY2oeT6uzvdm1j43sxOAImBFdMMTaRolDpHIugZ4yMw+AKzO8YeBNDNbDDwPXO3uFfW9gUi80Vh1ERFpFF1xiIhIoyhxiIhIoyhxiIhIoyhxiIhIoyhxiIhIoyhxiIhIoyhxiIhIoyhxiIhIo/x/rUDX3b1qF/wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"H1 error-estimate\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "\n",
    "plt.ion()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
